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The objective of the current investigation was to 
develop a simple turbulence model that will be applicable 
to propulsion flows having both bounded and unbounded 
regions by linking two existing algebraic turbulence 
models. The first is the Modified Mixing Length model, 
which is optimized for wall-bounded flows and has not 
previously been incorporated into the PARC code. The 
second is the Thomas model, the existing algebraic 
turbulence model in PARC, which has been used to 
calculate both bounded and unbounded turbulent flows, but 
was optimized for the latter. The following sections 
discuss both models, the method employed to link them 
into one model, and the validation of the resulting 
combination. 

Modified Mixing Length Model 
for Wall-Bounded Flows 

The Modified Mixing Length (MML) model is an 
algebraic turbulence model originally developed to analyze 
airflow over iced airfoils 3 with the ARC external flow 
Navier-Stokes code. 4 Potapczuk compared results obtained 
with this MML turbulence model and the Baldwin-Lomax 
turbulence model and found that the MML solutions more 
closely matched experimental data than did the Baldwin- 
Lomax solutions in separated flow situations. 

The MML model calculates turbulent viscosity for 
wall-bounded flows through Prandtl’s mixing length 
hypothesis 

A i, = P* 2 \(o\ (1) 

The mixing length / is defined as a function of distance 
from the wall 


and 
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The von Karman constant, K = 0.41, and van Driest 
constant, A + = 26, are the same as in the Baldwin-Lomax 
turbulence model, 5 and Cj and C 2 are empirical constants. 

The expressions for the length scale in a turbulent 
boundary layer given in Eqs. (2) and (3) were retained for 
the MML model in PARC. The terms C] and C 2 , 
however, were not restricted to being constants as in the 
original formulation; instead, their ratio was allowed to 
vary as a function of a local flow parameter, such as wall 
shear stress. An expression relating C\ and C 2 to the local 
wall shear stress r was constructed in the following 
manner. Equation (3) gives the capping value of the length 
scale /cap. which begins at the nondimensional position 
y + = Ci in the boundary layer and continues outward into 
the free stream. If this maximum length scale at a given 
position is related to the local boundary layer thickness 8 
through a constant B (typically B = 0.09), then the capping 
length scale is given as 


/ CAP = B 8 (6) 

Empirical correlations for 8 and r in Eqs. (7) and (8) are 
taken from Ref. 6 and are appropriate for a flat plate in the 
range of Reynolds numbers being considered; 
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Thus, an expression for the local capping length scale may 
( 4 ) be constructed as a function of the local skin friction or 
wall shear stress instead of boundary layer thickness: 


2 


(9) 


by firsl determining a length scale 


= (3.31'10' 6 )flc; 3 y 


Similar expressions that define the capping length scale 
as a function of boundary layer thickness, such as those of 
the Cebeci-Smith model, have been shown to be inadequate 
for separated flow situations. The edge of the boundary 
layer is difficult to determine in separated cases because of 
the recirculating flow. Equation (9) is defined 
independently of boundary layer thickness to improve the 
capability of PARC in predicting separated flows. For 
separated flows, Eq. (9) still will predict very large 
maximum length scales as the wall shear stress becomes 
small. To avoid this problem, Potapczuk^ used a weighted 
average for the shear stress 


t = .l|rj+ .2 |tJ+ .4|rJ+ .2 |t,J + .l|r,J a0 ) 


The subscripts in Eq. (10) refer to the computational grid 
point locations along the wall bounding the turbulent flow. 

Equations (3) and (9) may be used to give the 
desired relation among C i , C 2 , and shear stress (skin 
friction) 


— = (3.31 KT 6 /— ) C ' 35 01) 
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Equation (11) shows that if B and K are constants, the ratio 
of C\ to C 2 is a function of the shear stress. If either C\ 
or C 2 remains a constant, the other will be simply 
determined by Eq. (11) and will still allow their ratio to be 
a function of a local flow parameter such as shear stress. 

Thomas Model for Unbounded Flows 

The standard algebraic turbulence model in PARC 
is based on the work of Thomas . 7 Although the Thomas 
model in PARC can be applied to all flows, it was 
optimized for free shear layers, which are unbounded flows 
where, for example, a jet may be mixing with a slower 
flow. Its capability for calculating wall-bounded turbulent 
flows has been questioned. The unbounded part of the 
Thomas model calculates eddy viscosity in free shear layers 


t = e _o [MhD - MN)1 ( 12 ) 
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where uj is the velocity in a section, co c is the maximum 
vorticity in the section, and i 0 is a constant. Prandtfs 
mixing length hypothesis is used to determine turbulent 
viscosity as was done for the MML- bounded flow model. 
The Thomas-bounded flow model (not used for the current 
model) is detailed in Ref. 7. 

Linking of Bounded and Unbounded Models 

As mentioned previously, many propulsion flow 
cases will have both turbulent wall-bounded flow regions 
and unbounded regions. The current turbulence model will 
employ the MML model for the former and the Thomas 
model for the latter. The two models must also be linked 
appropriately to provide a smooth transition from the 
bounded to unbounded regions in the flow field. Figure 1 
depicts a flow situation having both a wall boundary layer 
and a free shear layer. In the current model, the MML 
model will calculate turbulent viscosity from the wall out 
to the nondimensional position y + = C 3 , which is roughly 
the edge of the wall boundary layer and can be defined by 
the local MML capping length scale multiplied by a 
constant. The Thomas model will be used in the 
unbounded region from the nondimensional position y+ = 
C 4 (also defined as a constant times the local MML 
capping length) out into the free stream. In the transition 
region between C 3 and C 4 , a linear function is used to 
determine the turbulent viscosity 


^mml( Q " y ) + /*Th {y ' Ci) (i3) 


The resulting combined model will be referred to as the 
MMLT model for the rest of this discussion. 


Validation of the MMLT Model 

Calculati on of Flow Over a Flat Plate 

The PARC code with the MMLT model was 
applied to a flow over a flat plate (M = 0.2) to determine 
appropriate values for C \ , C 2 , C 3 , and C 4 . Since C\ and 
C 2 are related by Eq. (11), five values of C 2 were 
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determined by selecting five different positions in the wall 
boundary layer (CyS+) where the MML capping length 
scale is reached. The parameter Q/S+ is the ratio of the 
distance from the wall where /cap is reached to the 
boundary layer thickness and 6+ is the nondimensional 
boundary layer thickness. Equations (3) and (6) can be used 
to obtain a relation between C 2 and the theoretical Cj/8 + . 



Note that the boundary layer thickness is used here only to 
optimize C 2 . As mentioned previously, the MMLT model 
does not use the boundary layer thickness to calculate 
turbulent length scales or turbulent viscosity. The ratio 
Cy 8+ was varied from 0.3 to 0.7 to determine both an 
optimal C 2 and the sensitivity of the model to this 
parameter. The beginning of the transition region, C 3 , was 
set as the nominal edge of the boundary layer, 10/ + CAP , 
and the end of the transition region, C 4 , was set to 
20/ + cap- 

Figures 2 and 3 show the effect that variation of 
Q/8+ had on the flat plate boundary layer predictions. The 
comparison of MMLT shear stress predictions to 
experimental data of Weighardt and Tillman 8 in Fig. 2 
indicates that the MMLT model matches the data best for 
larger values of C 1 /8 + . The MMLT boundary layer 
velocity profiles are compared to experimental data from 
Ref. 9 at a plate Reynolds number of 4 000 000 in Fig. 3a 
and at a plate Reynolds number of 10 000 000 in Fig. 3b. 
The MMLT solutions demonstrate no significant 
differences among the different values of Q/8+. Other 
combinations of C 3 and C 4 were examined and did not 
change the solutions provided 10/+ CAP ^ C 3 < C 4 . This 
was expected because C 3 and C 4 are the boundaries of the 
transition region from the wall-bounded MML model to the 
unbounded part of the Thomas model, and this flat plate 
flow has no free shear layer (unbounded) region. The 
following section discusses optimization of C 3 and C 4 for 
flow over a backward- facing step, a benchmark case for 
separating and reattaching flow. 

This same flat plate case was analyzed with the 
three turbulence models already available in PARC and the 
results were compared with the MMLT predictions using 
Cj/ 8+ = 0.7 (corresponding to C 2 = 3.18), C 3 = 10/ + C ap> 
and C 4 = 20/+ cap . The first of the three models was the 
Thomas model, the previously mentioned standard algebraic 
turbulence model in PARC. The second was the Baldwin- 
Lomax model recently installed by Sirbaugh. 10 The third 
and only two-equation model was Chien's low Reynolds 


number k-e model 11 modified by Nichols. 12 Figure 4 
shows a comparison of the shear stress predicted by these 
models with experimental data and Fig. 5 shows the 
boundary layer velocity profile comparisons at plate 
Reynolds numbers of 4 000 000 and 10 000 000 
respectively. Only the original Thomas model results 
disagree strongly with the data. 

The convergence histories for these flat plate cases 
are given in Table 1. The algebraic turbulence model 
solutions all took less than half the CPU time for the k-e 
solution, although all of the cases used the same 1 1 1 by 81 
grid. 

Calculation of Flow Over a Backward-Facing Step 

Flow over a backward-facing step has also been 
analyzed with the PARC code and the MMLT model in 
order to optimize C 3 and C 4 and to determine the model’s 
capability to calculate a separating and reattaching shear 
flow. Several MMLT backward-facing step cases were 
calculated (for varying Cj/8 + , C 3 and C 4 ) and compared. 
These PARC calculations modeled the experiment 
conducted by Driver and Seegmiller 13 for a flow over a 
backward-facing step with M = 0.128 upstream of the step. 
The experiment was conducted in a low-speed wind tunnel 
with an expansion ratio (tunnel height after step to tunnel 
height before step) of 1.125. The experimental geometry 
is shown in Fig. 6. The velocity profile measured at a 
position 4 step-heights upstream of the step in this 
experiment was used as the inflow boundary condition for 
the PARC calculations. 

The following four combinations of C 3 and C 4 
were used with Ci/5+ = 0.7 and Q/8+ = 0.3 (eight total 
combinations) for the MMLT cases: 

(1) C 3 = 10/+CAP* C 4 = 15/ + CAP 

(2) C 3 = 10/+ CA p, C 4 = 20/ + cap 

(3) C 3 = 10/ + c A p> C 4 = 40/ + c A p 

(4) C 3 = 20/ + c AP , C 4 = 40/ + cap 

Shear stress predictions for the Cj/8* = 0.7 cases 
appear in Fig. 7(a) and for the Ci/8 + = 0.3 cases in Fig. 
7(b). All skin friction results undershoot the experimental 
data in the separated region just downstream of the step, 
and all MMLT predictions of reattachment length (location 
where C f = 0) are larger than the experimental value. The 
PARC MMLT reattachment predictions are given in Table 
2. Eaton and Johnston 14 report that negative skin friction 
coefficients as large as -0.0012 are common for backward- 
facing step flows. This is still smaller than the MMLT 
predictions but larger than the experimental data of Driver 
and Seegmiller. The most accurate reattachment 
predictions are provided by the solution with C 3 = 10/ + CAP 
and either C 4 = 15/ + c AP or C 4 = 20/ + (; AP . The MMLT 
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skin friction predictions in the reattached region 
downstream of the step match the experimental data well. 
There is much less variation among the various 
combinations of C 3 and C 4 for cases with cy5+= 0.7 than 
with Q/8 + = 0.3. 

Velocity profiles at X/H = 3, 6, 12, and 20 appear 
in Figs. 8 and 9 for cases with Q/8+ = 0.7 and Cy8 + = 
0.3, respectively. The MMLT profiles overestimate the 
magnitude of the reversed velocity in the separated region 
(X/H = 3 and X/H = 6) with the largest discrepancy for C 4 
= 40/+ CAP . This corresponds to the skin friction results. 
The experimental velocities at X/H = 6 indicate that the 
flow has already reattached at this position, although the 
experimental skin friction results indicate the flow is still 
separated, as Driver and Seegmiller report. 

As was done for the flat plate case, backward- 
facing step calculations were also obtained with PARCs 
three other turbulence models for comparison with the 
MMLT model. The MMLT solution with Q/8+ = 0.7, C 3 
- 10/ + CAP , and C 4 = 20 £ + cap was used for comparison 
with the Thomas, Baldwin-Lomax, and k-e solutions. 
Figure 10 shows a comparison of the skin friction 
predictions for these models and Table 3 gives the 
reattachment length predictions. Compared with the 
experimental data, the Thomas solution shows a much 
smaller negative skin friction before the reattachment and a 
much smaller positive skin friction after. The location of 
the reattachment is also severely overestimated. 

The Baldwin-Lomax solution demonstrates the 
largest undershoot of the skin friction before the 
reattachment and largest overshoot al ter the reattachment of 
any of die models. The Baldwin-Lomax model behaves 
poorly for separated flows because the key function in the 
model, which describes the product of vorticity and length 
scale, is not easily determined for such situations. The 
Baldwin-Lomax model also predicts a smaller reattachment 
length than the experimental data. 

It was anticipated that die k-e solution would 
provide the best match to the experimental data. The k-e 
solution does predict nearly the same reattachment length 
as the experimental data, but significantly overestimates 
die magnitude of the skin friction before and after the 
reattachment. The k-e skin friction prediction is the only 
one to reach a peak after the reattachment and then become 
smaller farther downstream. Avva, Smith, and Singhal 15 
report the same skin friction behavior when applying a low 
Reynolds number k-e model to the same backward-facing 
step of Driver and Seegmiller. They found that die k- 
e solution is very sensitive to the grid packing in die inner 
layer (y + < 30) by varying the number of grid points in 


this region from 5 to 30. The current study used the same 
grid for all the PARC backward- facing step cases with 
approximately 18 grid points in the inner layer. 

As mentioned previously, the MMLT model 
undershoots the skin friction before the reattachment and 
predicts a larger reattachment length than was 
experimentally determined. However, the MMLT soludon 
is the only one to closely match the experimental skin 
friction data downstream of the reattachment position. 

Figure 11 presents a comparison of the velocity 
profiles at X/H = 3, 6, 12, and 20. The Thomas solution 
overall shows the poorest agreement with the experimental 
data as was the case with the flat plate examination. The 
other three models match the experimental data more 
accurately, but there is a large variation in their velocity 
profiles near the wall (Y/H < 1). Away from the wall, the 
Baldwin-Lomax solution predicts a larger free-stream 
velocity than do the other models or the experimental data. 

The convergence histories for the backward-facing 
step cases are given in Table 4. These cases took much 
longer to converge than the flat plate cases. This was most 
likely due to the following: (1) the increased complexity of 
the separating and reattaching flow that caused the 
maximum allowable time step in PARC to be 30 times 
smaller than for the flat plate and (2) the PARC code’s 
convergence rate becoming very slow for flows with a free- 
stream Mach number near 0.1 or smaller (0.128 for the 
backward-facing step). The k-e solution took more than 
twice the iterations and more than seven times the CPU 
time to come to convergence than any of die algebraic 
model solutions. 

Conclusions 

The flat plate and backward-facing step flow 
computations with the PARC code have been valuable 
both in developing the MMLT model from its two existing 
algebraic turbulence model components and in assessing its 
capabilities and those of the other turbulence models in 
PARC. 

All the turbulence models, except PARC’s 
standard algebraic turbulence model, the Thomas model, 
provided accurate skin friction and boundary layer velocity 
profile predictions for die flat plate flow. None of the 
models agreed very well with the experimental skin friction 
and velocity profile data for the backward-facing step case 
in the separated region. Downstream of the reattachment, 
all the turbulence models, except for Thomas, show fair 
agreement with the experimental velocity profiles. The 
MMLT skin friction results match the experimental data 
downstream of the reattachment much better than the odier 
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models, including the k-e model which took seven times 
more CPU time than any of the algebraic models to 
converge. 

Other flow cases will be investigated to determine 
the PARC code’s capability to provide accurate propulsion 
flow predictions with the MMLT model. A single flow 
plug nozzle is being constructed and will be tested by 
NASA Langley Research Center to provide extensive data 
for code validation. This nozzle has been investigated with 
PARC using the Thomas and k-e models and will be 
investigated with PARC using the MMLT model. 
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Table 1. Convergence histories for flat plate solutions. 


Table 3. Comparison of reattachment positions of 
PARC's other turbulence model solutions. 


MODEL 

ITERATIONS 

CRAY Y/MP 
CPU TIME (s) 

MMLT 

4000 

500 

Thomas 

4000 

450 

Baldwin-Lomax 

4000 

500 

Chien k-e 

5000 

1100 


Table 2. Comparison of reattachment positions of the 
MMLT solutions. 


CASE 

REATTACHMENT 
POSITION 
(STEP HEIGHTS, H) 

Driver-Seegmiller Data 

6.250 

(C1/8+ = .7, C3 = 10/ + QSJ>> 
C4 = 15/ + cap) 

7.353 

(C]/S+ = .7, C3 = 10/-CAP, 
C4 = 20 / -cap) 

7.416 

(Cj/8+ = .7, C 3 = 10/ -cap. 
C4 = 40/ -cap) 

7.566 

(C|/8+ = .7, C 3 = 20/ -cap. 
C4 = 40/+cap) 

7.749 

(C[/8+ — .7, C 3 — 10/ + cap» 
C4 = 15/ -cap) 

7.143 

(C,/8+ = .7, C 3 = 10/ -cap, 
C 4 = 20/ -cap) 

7.429 

(Cj/8+ = .7, C3 = 10/-CAP. 
C4 = 40/ -cap) 

7.746 


CASE 

REATTACHMENT 
POSITION 
(STEP HEIGHTS, H) 

Driver-Seegmiller Data 

6.250 

MMLT (Q/8+ = .7, 

C3 = IO/+CAP.C4 = 20/^cap) 

7.416 

Thomas 

12.281 

Baldwin-Lomax 

5.410 

Chien k-e 

6.256 


Table 4. Convergence histories for backward- facing step 
solutions. 


MODEL 

ITERATIONS 

CRAY Y/MP 
CPU TIME (s) 

MMLT 

45,000 

7100 

Thomas 

40,000 

6000 

Baldwin-Lomax 

45,000 

7000 

Chien k-e 

110,000 

50,000 


(C1/8+ = .7, C3 = 20/ + cap, 
C 4 = 40£+c AP ) 


8.786 


















Local skin friction, 


t 

THOMAS UNBOUNDED REGION 




2000 4000 6000 8000 10000 12000 

Reynolds number base on momentum 
thickness, Re 0 


Figure 2. — MMLT flat plate skin friction predictions 
for varying C-|/8 + . 



Position in boundary layer, y/S 


(b) Re x = 1 0 7 . 

Figure 3. — MMLT flat plate boundary layer velocity 
profiles for varying C 1 /5 + . 



Reynolds number based on momentum 
thickness, Re 0 

Figure 4. — Flat plate skin friction predictions using 
the different turbulence models in PARC. 
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Position in boundary layer, y/8 
(b) Re x = 10 7 . 


MMLT 

c 3 = 1 0/ + CA p, C 4 = 1 5 / + cap 

C 3 = 1 0/ + CA p, C 4 = 20/ + cap 
C 3 = lOZ+^p, C 4 = 40 / + cap 
C 3 = 20/ + cap , C 4 = 40 / + cap 
° Driver-Seegmiller 



c 


Figure 5. — Flat plate boundary layer velocity 
profiles using the different turbulence models 
in PARC. 




Axial position, X/H 
(b) C 1 /8 + = .3. 

Figure 7. — MMLT backward -facing step skin 
friction predictions with varying C 3 and C 4 . 


Figure 6. — Geometry of backward-facing step in experiment 
of Driver and Seegmiller. 
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Vertical position, Y/H 


MMLT (C 3 = 1 0r*' CAPf C 4 = 20/+ CAP ) 

Thomas 

Baldwin- Lomax 
K-Epsilon 

o Driver-Seegmiller 



Axial position, X/H 


Figure 1 0. — Backward-facing step skin friction 
predictions using the different turbulence 
models in PARC. 


MMLT (C 3 = 1 0/ + CAPt C 4 = 20/ + CAP ) 

Thomas 

Baldwin-Lomax 

K-Epsilon 

o Driver-Seemilier 






(b) X/H = 6. (d) X/H * 20. 


u/u ref 

Figure 1 1 .-Backward-facing step velocity profile predictions using the different turbulence 
models in PARC. 
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